1 Introduction and summary

2 Loading the data

train <- read.csv("train.csv"); test <- read.csv("test.csv")

3 Summarizing the data

summary(train)
##        Id           MSSubClass       MSZoning     LotFrontage    
##  Min.   :   1.0   Min.   : 20.0   C (all):  10   Min.   : 21.00  
##  1st Qu.: 365.8   1st Qu.: 20.0   FV     :  65   1st Qu.: 59.00  
##  Median : 730.5   Median : 50.0   RH     :  16   Median : 69.00  
##  Mean   : 730.5   Mean   : 56.9   RL     :1151   Mean   : 70.05  
##  3rd Qu.:1095.2   3rd Qu.: 70.0   RM     : 218   3rd Qu.: 80.00  
##  Max.   :1460.0   Max.   :190.0                  Max.   :313.00  
##                                                  NA's   :259     
##     LotArea        Street      Alley      LotShape  LandContour
##  Min.   :  1300   Grvl:   6   Grvl:  50   IR1:484   Bnk:  63   
##  1st Qu.:  7554   Pave:1454   Pave:  41   IR2: 41   HLS:  50   
##  Median :  9478               NA's:1369   IR3: 10   Low:  36   
##  Mean   : 10517                           Reg:925   Lvl:1311   
##  3rd Qu.: 11602                                                
##  Max.   :215245                                                
##                                                                
##   Utilities      LotConfig    LandSlope   Neighborhood   Condition1  
##  AllPub:1459   Corner : 263   Gtl:1382   NAmes  :225   Norm   :1260  
##  NoSeWa:   1   CulDSac:  94   Mod:  65   CollgCr:150   Feedr  :  81  
##                FR2    :  47   Sev:  13   OldTown:113   Artery :  48  
##                FR3    :   4              Edwards:100   RRAn   :  26  
##                Inside :1052              Somerst: 86   PosN   :  19  
##                                          Gilbert: 79   RRAe   :  11  
##                                          (Other):707   (Other):  15  
##    Condition2     BldgType      HouseStyle   OverallQual    
##  Norm   :1445   1Fam  :1220   1Story :726   Min.   : 1.000  
##  Feedr  :   6   2fmCon:  31   2Story :445   1st Qu.: 5.000  
##  Artery :   2   Duplex:  52   1.5Fin :154   Median : 6.000  
##  PosN   :   2   Twnhs :  43   SLvl   : 65   Mean   : 6.099  
##  RRNn   :   2   TwnhsE: 114   SFoyer : 37   3rd Qu.: 7.000  
##  PosA   :   1                 1.5Unf : 14   Max.   :10.000  
##  (Other):   2                 (Other): 19                   
##   OverallCond      YearBuilt     YearRemodAdd    RoofStyle   
##  Min.   :1.000   Min.   :1872   Min.   :1950   Flat   :  13  
##  1st Qu.:5.000   1st Qu.:1954   1st Qu.:1967   Gable  :1141  
##  Median :5.000   Median :1973   Median :1994   Gambrel:  11  
##  Mean   :5.575   Mean   :1971   Mean   :1985   Hip    : 286  
##  3rd Qu.:6.000   3rd Qu.:2000   3rd Qu.:2004   Mansard:   7  
##  Max.   :9.000   Max.   :2010   Max.   :2010   Shed   :   2  
##                                                              
##     RoofMatl     Exterior1st   Exterior2nd    MasVnrType    MasVnrArea    
##  CompShg:1434   VinylSd:515   VinylSd:504   BrkCmn : 15   Min.   :   0.0  
##  Tar&Grv:  11   HdBoard:222   MetalSd:214   BrkFace:445   1st Qu.:   0.0  
##  WdShngl:   6   MetalSd:220   HdBoard:207   None   :864   Median :   0.0  
##  WdShake:   5   Wd Sdng:206   Wd Sdng:197   Stone  :128   Mean   : 103.7  
##  ClyTile:   1   Plywood:108   Plywood:142   NA's   :  8   3rd Qu.: 166.0  
##  Membran:   1   CemntBd: 61   CmentBd: 60                 Max.   :1600.0  
##  (Other):   2   (Other):128   (Other):136                 NA's   :8       
##  ExterQual ExterCond  Foundation  BsmtQual   BsmtCond    BsmtExposure
##  Ex: 52    Ex:   3   BrkTil:146   Ex  :121   Fa  :  45   Av  :221    
##  Fa: 14    Fa:  28   CBlock:634   Fa  : 35   Gd  :  65   Gd  :134    
##  Gd:488    Gd: 146   PConc :647   Gd  :618   Po  :   2   Mn  :114    
##  TA:906    Po:   1   Slab  : 24   TA  :649   TA  :1311   No  :953    
##            TA:1282   Stone :  6   NA's: 37   NA's:  37   NA's: 38    
##                      Wood  :  3                                      
##                                                                      
##  BsmtFinType1   BsmtFinSF1     BsmtFinType2   BsmtFinSF2     
##  ALQ :220     Min.   :   0.0   ALQ :  19    Min.   :   0.00  
##  BLQ :148     1st Qu.:   0.0   BLQ :  33    1st Qu.:   0.00  
##  GLQ :418     Median : 383.5   GLQ :  14    Median :   0.00  
##  LwQ : 74     Mean   : 443.6   LwQ :  46    Mean   :  46.55  
##  Rec :133     3rd Qu.: 712.2   Rec :  54    3rd Qu.:   0.00  
##  Unf :430     Max.   :5644.0   Unf :1256    Max.   :1474.00  
##  NA's: 37                      NA's:  38                     
##    BsmtUnfSF       TotalBsmtSF      Heating     HeatingQC CentralAir
##  Min.   :   0.0   Min.   :   0.0   Floor:   1   Ex:741    N:  95    
##  1st Qu.: 223.0   1st Qu.: 795.8   GasA :1428   Fa: 49    Y:1365    
##  Median : 477.5   Median : 991.5   GasW :  18   Gd:241              
##  Mean   : 567.2   Mean   :1057.4   Grav :   7   Po:  1              
##  3rd Qu.: 808.0   3rd Qu.:1298.2   OthW :   2   TA:428              
##  Max.   :2336.0   Max.   :6110.0   Wall :   4                       
##                                                                     
##  Electrical     X1stFlrSF      X2ndFlrSF     LowQualFinSF    
##  FuseA:  94   Min.   : 334   Min.   :   0   Min.   :  0.000  
##  FuseF:  27   1st Qu.: 882   1st Qu.:   0   1st Qu.:  0.000  
##  FuseP:   3   Median :1087   Median :   0   Median :  0.000  
##  Mix  :   1   Mean   :1163   Mean   : 347   Mean   :  5.845  
##  SBrkr:1334   3rd Qu.:1391   3rd Qu.: 728   3rd Qu.:  0.000  
##  NA's :   1   Max.   :4692   Max.   :2065   Max.   :572.000  
##                                                              
##    GrLivArea     BsmtFullBath     BsmtHalfBath        FullBath    
##  Min.   : 334   Min.   :0.0000   Min.   :0.00000   Min.   :0.000  
##  1st Qu.:1130   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:1.000  
##  Median :1464   Median :0.0000   Median :0.00000   Median :2.000  
##  Mean   :1515   Mean   :0.4253   Mean   :0.05753   Mean   :1.565  
##  3rd Qu.:1777   3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:2.000  
##  Max.   :5642   Max.   :3.0000   Max.   :2.00000   Max.   :3.000  
##                                                                   
##     HalfBath       BedroomAbvGr    KitchenAbvGr   KitchenQual
##  Min.   :0.0000   Min.   :0.000   Min.   :0.000   Ex:100     
##  1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:1.000   Fa: 39     
##  Median :0.0000   Median :3.000   Median :1.000   Gd:586     
##  Mean   :0.3829   Mean   :2.866   Mean   :1.047   TA:735     
##  3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:1.000              
##  Max.   :2.0000   Max.   :8.000   Max.   :3.000              
##                                                              
##   TotRmsAbvGrd    Functional    Fireplaces    FireplaceQu   GarageType 
##  Min.   : 2.000   Maj1:  14   Min.   :0.000   Ex  : 24    2Types :  6  
##  1st Qu.: 5.000   Maj2:   5   1st Qu.:0.000   Fa  : 33    Attchd :870  
##  Median : 6.000   Min1:  31   Median :1.000   Gd  :380    Basment: 19  
##  Mean   : 6.518   Min2:  34   Mean   :0.613   Po  : 20    BuiltIn: 88  
##  3rd Qu.: 7.000   Mod :  15   3rd Qu.:1.000   TA  :313    CarPort:  9  
##  Max.   :14.000   Sev :   1   Max.   :3.000   NA's:690    Detchd :387  
##                   Typ :1360                               NA's   : 81  
##   GarageYrBlt   GarageFinish   GarageCars      GarageArea     GarageQual 
##  Min.   :1900   Fin :352     Min.   :0.000   Min.   :   0.0   Ex  :   3  
##  1st Qu.:1961   RFn :422     1st Qu.:1.000   1st Qu.: 334.5   Fa  :  48  
##  Median :1980   Unf :605     Median :2.000   Median : 480.0   Gd  :  14  
##  Mean   :1979   NA's: 81     Mean   :1.767   Mean   : 473.0   Po  :   3  
##  3rd Qu.:2002                3rd Qu.:2.000   3rd Qu.: 576.0   TA  :1311  
##  Max.   :2010                Max.   :4.000   Max.   :1418.0   NA's:  81  
##  NA's   :81                                                              
##  GarageCond  PavedDrive   WoodDeckSF      OpenPorchSF     EnclosedPorch   
##  Ex  :   2   N:  90     Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
##  Fa  :  35   P:  30     1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:  0.00  
##  Gd  :   9   Y:1340     Median :  0.00   Median : 25.00   Median :  0.00  
##  Po  :   7              Mean   : 94.24   Mean   : 46.66   Mean   : 21.95  
##  TA  :1326              3rd Qu.:168.00   3rd Qu.: 68.00   3rd Qu.:  0.00  
##  NA's:  81              Max.   :857.00   Max.   :547.00   Max.   :552.00  
##                                                                           
##    X3SsnPorch      ScreenPorch        PoolArea        PoolQC    
##  Min.   :  0.00   Min.   :  0.00   Min.   :  0.000   Ex  :   2  
##  1st Qu.:  0.00   1st Qu.:  0.00   1st Qu.:  0.000   Fa  :   2  
##  Median :  0.00   Median :  0.00   Median :  0.000   Gd  :   3  
##  Mean   :  3.41   Mean   : 15.06   Mean   :  2.759   NA's:1453  
##  3rd Qu.:  0.00   3rd Qu.:  0.00   3rd Qu.:  0.000              
##  Max.   :508.00   Max.   :480.00   Max.   :738.000              
##                                                                 
##    Fence      MiscFeature    MiscVal             MoSold      
##  GdPrv:  59   Gar2:   2   Min.   :    0.00   Min.   : 1.000  
##  GdWo :  54   Othr:   2   1st Qu.:    0.00   1st Qu.: 5.000  
##  MnPrv: 157   Shed:  49   Median :    0.00   Median : 6.000  
##  MnWw :  11   TenC:   1   Mean   :   43.49   Mean   : 6.322  
##  NA's :1179   NA's:1406   3rd Qu.:    0.00   3rd Qu.: 8.000  
##                           Max.   :15500.00   Max.   :12.000  
##                                                              
##      YrSold        SaleType    SaleCondition    SalePrice     
##  Min.   :2006   WD     :1267   Abnorml: 101   Min.   : 34900  
##  1st Qu.:2007   New    : 122   AdjLand:   4   1st Qu.:129975  
##  Median :2008   COD    :  43   Alloca :  12   Median :163000  
##  Mean   :2008   ConLD  :   9   Family :  20   Mean   :180921  
##  3rd Qu.:2009   ConLI  :   5   Normal :1198   3rd Qu.:214000  
##  Max.   :2010   ConLw  :   5   Partial: 125   Max.   :755000  
##                 (Other):   9

The outcome variable SalePrice is a complete variable, but there are a number of predictors with many missing values. Median house price is $163,000.

par(mfrow = c(1,2))
boxplot(train$SalePrice,main = "SalePrice")
boxplot(log10(train$SalePrice), main = "Log10(SalePrice)")

3.1 Understanding the impact of missing data

Which features have more than 30% missing data?

na.status <- is.na(train)
na.sum <- apply(na.status,2,sum)
names(na.sum) <- colnames(train)
mostly_missing <- which(na.sum > (0.3 * nrow(train)))
na.sum[mostly_missing]
##       Alley FireplaceQu      PoolQC       Fence MiscFeature 
##        1369         690        1453        1179        1406

5 variables have substantial amount of missing data. Do they have impact on the house prices?

3.1.1 Alley

Based on the data description: Type of alley access to property

   Grvl     Gravel
   Pave     Paved
   NA       No alley access
   
library(ggplot2)

ggplot(data = train, aes(x = Alley, y = log(SalePrice), fill= Alley))+
        geom_boxplot()

Alley might be important. Gravel homes have significantly lower price. Convert NAs to “NoAlley” to make more sense.

train$Alley <- as.character(train$Alley)
train$Alley[which(is.na(train$Alley))] <- "NoAlley"
train$Alley <- factor(train$Alley)

#Also transform the test set
test$Alley <- as.character(test$Alley)
test$Alley[which(is.na(test$Alley))] <- "NoAlley"
test$Alley <- factor(test$Alley)

3.1.2 FireplaceQu

Based on the data description:

FireplaceQu: Fireplace quality

   Ex       Excellent - Exceptional Masonry Fireplace
   Gd       Good - Masonry Fireplace in main level
   TA       Average - Prefabricated Fireplace in main living area or Masonry Fireplace in basement
   Fa       Fair - Prefabricated Fireplace in basement
   Po       Poor - Ben Franklin Stove
   NA       No Fireplace
ggplot(data = train, aes(x = FireplaceQu, y = log10(SalePrice), fill= FireplaceQu))+
        geom_boxplot()

Let’s look at through a simple linear model:

par(mfrow=c(1,2))
hist(train$SalePrice)
hist(log10(train$SalePrice)) # Looks more gaussian

FirePlaceFit <- lm(log10(train$SalePrice) ~ FireplaceQu, data = train)
summary(FirePlaceFit)
## 
## Call:
## lm(formula = log10(train$SalePrice) ~ FireplaceQu, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38691 -0.09351 -0.01374  0.09604  0.57966 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.50252    0.03013 182.644  < 2e-16 ***
## FireplaceQuFa -0.28729    0.03959  -7.256 9.82e-13 ***
## FireplaceQuGd -0.17957    0.03106  -5.781 1.08e-08 ***
## FireplaceQuPo -0.40442    0.04469  -9.050  < 2e-16 ***
## FireplaceQuTA -0.21003    0.03126  -6.719 3.58e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1476 on 765 degrees of freedom
##   (690 observations deleted due to missingness)
## Multiple R-squared:  0.118,  Adjusted R-squared:  0.1134 
## F-statistic: 25.59 on 4 and 765 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(FirePlaceFit)

It appears that all levels of this feature also have significantly different mean Sales prices. Again, converting NAs to NoFireplace to be more informative.

train$FireplaceQu <- as.character(train$FireplaceQu)
train$FireplaceQu[which(is.na(train$FireplaceQu))] <- "NoFireplace"
train$FireplaceQu <- factor(train$FireplaceQu)

#Also transform the test set
test$FireplaceQu <- as.character(test$FireplaceQu)
test$FireplaceQu[which(is.na(test$FireplaceQu))] <- "NoFireplace"
test$FireplaceQu <- factor(test$FireplaceQu)

3.1.3 PoolQC

Based on the data description:

    PoolQC: Pool quality
           Ex       Excellent
           Gd       Good
           TA       Average/Typical
           Fa       Fair
           NA       No Pool
   
library(ggplot2)
ggplot(data = train, aes(x = PoolQC, y = log10(SalePrice), fill = PoolQC))+
        geom_boxplot()

Note that actually very few houses have pools. Houses with excellent pool condition have significantly higher sales prices. We will also keep this feature, converting NAs to “NoPool”:

train$PoolQC <- as.character(train$PoolQC)
train$PoolQC[which(is.na(train$PoolQC))] <- "NoPool"
train$PoolQC <- factor(train$PoolQC)

#Also transform the test set
test$PoolQC <- as.character(test$PoolQC)
test$PoolQC[which(is.na(test$PoolQC))] <- "NoPool"
test$PoolQC <- factor(test$PoolQC)

3.1.4 Fence

Based on the data description:

    Fence: Fence quality
           GdPrv    Good Privacy
           MnPrv    Minimum Privacy
           GdWo     Good Wood
           MnWw     Minimum Wood/Wire
           NA       No Fence 
library(ggplot2)
ggplot(data = train, aes(x = Fence, y = log10(SalePrice), fill = Fence))+
        geom_boxplot()

FenceFit <- lm(log10(train$SalePrice) ~ Fence, data = train)
summary(FenceFit)
## 
## Call:
## lm(formula = log10(train$SalePrice) ~ Fence, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57500 -0.06747 -0.00520  0.05611  0.72551 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.23634    0.01816 288.272  < 2e-16 ***
## FenceGdWo   -0.11851    0.02628  -4.510 9.58e-06 ***
## FenceMnPrv  -0.08969    0.02131  -4.210 3.46e-05 ***
## FenceMnWw   -0.11305    0.04582  -2.467   0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1395 on 277 degrees of freedom
##   (1179 observations deleted due to missingness)
## Multiple R-squared:  0.08211,    Adjusted R-squared:  0.07217 
## F-statistic:  8.26 on 3 and 277 DF,  p-value: 2.781e-05

Interestingly, this feature might have a complex relationship, if any at all, with the response variable. The effects of each of the level deemed as significant, so we will keep this feature at this point as well.

train$Fence <- as.character(train$Fence)
train$Fence[which(is.na(train$Fence))] <- "NoFence"
train$Fence <- factor(train$Fence)

#Also transform the test set
test$Fence <- as.character(test$Fence)
test$Fence[which(is.na(test$Fence))] <- "NoFence"
test$Fence <- factor(test$Fence)

3.1.5 MiscFeature

Based on the data description:

    MiscFeature: Miscellaneous feature not covered in other categories
           Elev     Elevator
           Gar2     2nd Garage (if not described in garage section)
           Othr     Other
           Shed     Shed (over 100 SF)
           TenC     Tennis Court
           NA       None
           

We will deal with the remaining missing value containing features in the next stages of our analysis.

table(train$MiscFeature)
## 
## Gar2 Othr Shed TenC 
##    2    2   49    1
ggplot(data = train, aes(x = MiscFeature, y = log10(SalePrice), fill = MiscFeature))+
        geom_boxplot()

Again, very few houses have this features, and not sure if these qualities have a real impact on the house value amongst other measured features. Neverthless, I will maintain this feature as well. Just converting NAs to “NoMiscFeature” to be more clear:

train$MiscFeature <- as.character(train$MiscFeature)
train$MiscFeature[which(is.na(train$MiscFeature))] <- "NoMiscFeature"
train$MiscFeature <- factor(train$MiscFeature)

#Also transform the test set
test$MiscFeature <- as.character(test$MiscFeature)
test$MiscFeature[which(is.na(test$MiscFeature))] <- "NoMiscFeature"
test$MiscFeature <- factor(test$MiscFeature)

4 Developing expectations from data

Here I will explore the data further, in order to develope a better understanding about the relationships between features and the home prices.

4.1 Continuous features: watching for strong predictors and potential collinearities

Let’s look at the continuous features in the data set:

cont.var <-NULL
for(i in 1:ncol(train)){
        if(class(train[,i]) == "integer" | class(train[,i]) == "numeric"){
                cont.var <- c(cont.var,i)
        }
}

4.1.1 Part1: transforming LotArea to make it more gaussian

pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[1:10]], cex = 0.05, pch = 19, col = "navy")

  • We need to remove ID, since this is just an identifier for the houses, has no predictive value as expected.
  • OverallQual and OverallCond look strong predictors, but might have some collinearity.
  • Yearbuilt and YearRemodAdd also look good predictors but most likely contain redundant information.
  • For the LOT and SF features, it could be better to look at after log transformation, since the data looks skewed.
par(mfrow=c(3,2))
hist(train$LotArea, breaks = 50, col = "navy");hist(log10(train$LotArea+1),breaks=50,col = "navy")
hist(train$MasVnrArea, breaks = 50,col = "navy");hist(log10(train$MasVnrArea+1),breaks=50,col = "navy")
hist(train$BsmtFinSF1, breaks = 50,col = "navy");hist(log10(train$BsmtFinSF1+1),breaks=50,col = "navy")

Indeed one of these features (LotArea) look more gaussian when it is log10 transformed.

pairs(log10(train$SalePrice) ~ log10(train$LotArea+1) + log10(train$MasVnrArea+1) + log10(train$BsmtFinSF1+1),cex = 0.1, pch = 19, col = "navy")

When transformed, lotArea seems to do a better job in explaining the variability in the SalePrice, other 2 variables however, have substantial 0 values and it is questionable how much this transformation might help.

Therefore, I will only transform the LotArea at this stage and leave the other features as they are, for feature selection steps.

4.1.3 Part3: Size of the garage and number of rooms are good predictors

pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[22:31]], cex = 0.05, pch = 19, col = "navy")

Arguably, these features are all related to the overall area of the house.

4.1.4 Part4: fairly weak predictors

pairs(log10(train$SalePrice) ~ ., data = train[,cont.var[32:38]], cex = 0.05, pch = 19, col = "navy")

Processing the identified features:

# Remove ID variable

train.ID <- train$Id #Hold it in a new variable
train <- train[, -1]

library(dplyr)

train <- dplyr::mutate(train, log10.LotArea = log10(LotArea))
train <- dplyr::select(train, -LotArea)

# Also transform in the test set
test <- dplyr::mutate(test, log10.LotArea = log10(LotArea))
test <- dplyr::select(test, -LotArea)

4.2 Categorical predictors: feature engineering to create dummy variables

Now let’s look at the categorical variables to see if certain ones are more powerful than others:

cat.var <-NULL
for(i in 1:ncol(train)){
        if(class(train[,i]) == "factor"){
                cat.var <- c(cat.var,i)
        }
}
length(cat.var)
## [1] 43

We have 43 categorical variables, some of them have only few levels. This is a challenge in many ways. It is not feasible to look each of them individually, therefore, we will first convert them to binary variables and try to use regularization or decision tree-based approaches to select and use most useful features.

# Write a function that gets a data set and converts all factor variables to dummy variables
convert.to.dummy <- function(data.set){
        cat.var <-NULL
        temp.data <- data.frame(1:nrow(data.set))
           for(i in 1:ncol(data.set)){
                if(class(data.set[,i]) == "factor"){
                        cat.var <- c(cat.var,i)
                        factor.levels <- levels(data.set[,i]) # Try to find a way to classify NA's as "NO" otherwise they generate problem downstream
                                # First check if there is any 'NA-level'
                                if(any(is.na(data.set[,i]))){
                                        dummy.vector = ifelse(is.na(data.set[,i]),1,0)
                                        dummy.vector <- data.frame(dummy.vector)
                                        colnames(dummy.vector)[1] = paste("NO",names((data.set)[i]),sep = ".")
                                        temp.data <- cbind(temp.data,dummy.vector)
                                }
                        
                                for(j in seq_along(factor.levels)){ # Then deal with normal factor levels
                                dummy.vector = ifelse(data.set[,i] == factor.levels[j],1,0)
                                
                                #Since we already dealt with NAs above
                                if(any(is.na(dummy.vector))){dummy.vector[is.na(dummy.vector)] <- 0} 
                                
                                dummy.vector <- data.frame(dummy.vector)
                                colnames(dummy.vector)[1] = paste(names((data.set)[i]),
                                                                  factor.levels[j],sep = ".")
                                temp.data <- cbind(temp.data,dummy.vector)
                                }
                }
           }
           #Remove the original categorical variables from data.set
           data.set <- data.set[,-cat.var]     
           #Add the dummy.variable set
           temp.data <- temp.data[,-1] # remove the unnecessary column
           data.set <- cbind(data.set,temp.data)
           
           return(data.set)     
}

# Keep test.IDs aside
test.ID <- test$Id

# Process the training set
training.processed <- convert.to.dummy(train)

# Process the test set
test.processed <- convert.to.dummy(test)

# Note that not all the levels are present in the test set, therefore we need to make the features same.
training.processed.SalePrice <- training.processed$SalePrice #Keep the original outcome variable 

consensus.features1 <- which(colnames(training.processed) %in% colnames(test.processed))
training.processed <- training.processed[,consensus.features1] #Keep the same features as in the test set

consensus.features2 <- which(colnames(test.processed) %in% colnames(training.processed))
test.processed <- test.processed[,consensus.features2] #Keep the same features as in the training set

identical(colnames(test.processed), colnames(training.processed))
## [1] TRUE
# Same features in both sets!

training.processed <- data.frame(Log.SalePrice = log(training.processed.SalePrice),training.processed) # Add the Log of the response variable
test.processed <- data.frame(ID = test.ID,test.processed) # Add the ID variable

4.3 Missing value imputation

We still have missing values in some of our features, which will be a problem down the road.

apply(is.na(training.processed), 2, sum)[apply(is.na(training.processed), 2, sum) != 0]
## LotFrontage  MasVnrArea GarageYrBlt 
##         259           8          81
apply(is.na(test.processed), 2, sum)[apply(is.na(test.processed), 2, sum) != 0]
##  LotFrontage   MasVnrArea   BsmtFinSF1   BsmtFinSF2    BsmtUnfSF 
##          227           15            1            1            1 
##  TotalBsmtSF BsmtFullBath BsmtHalfBath  GarageYrBlt   GarageCars 
##            1            2            2           78            1 
##   GarageArea 
##            1

4.4 GarageYrBlt

We have 81 missing cases in our training set, these are the houses actually they don’t have a garage.

plot(y =training.processed$Log.SalePrice, x =training.processed$GarageYrBlt, pch = 19, col = "navy", cex = 0.5)

It looks like this feature has some relationship with the outcome, as also we can generally expect that newer houses are going to be more expensive. Therefore, we should make an attempt to impute the values missing from this feature.

plot(training.processed$YearBuilt, training.processed$GarageYrBlt, pch = 19, col = "navy", cex = 0.5)

As we expected, for most of the houses the year that garage was built is the same as the year the house was built. Therefore, for imputation purpose it is sensible to use the YearBuilt feature to fill these issing values in the GarageYrBuilt feature. While this is not consistent with the fact that these houses do not contain any garage, since we captured that in the categorical features we engineered, at this point imputing the years seems to be more appropriate than removing these data points alltogether:

#Impute in the training set
training.processed$GarageYrBlt[is.na(training.processed$GarageYrBlt)] <- training.processed$YearBuilt[is.na(training.processed$GarageYrBlt)]


#Impute in the test set
test.processed$GarageYrBlt[is.na(test.processed$GarageYrBlt)] <- test.processed$YearBuilt[is.na(test.processed$GarageYrBlt)]

4.5 LotFrontage

Let’s look at this feature more closely:

plot(x = training.processed$LotFrontage, y = training.processed$Log.SalePrice, cex = 0.5, col = "navy" , pch = 19)
abline(lm(training.processed$Log.SalePrice ~ training.processed$LotFrontage))

Again this feature seems to have some degree of relationship with the response. We suspect that this feature has relationship with other features that relate to the Lot or the general area of the house:

cor(training.processed$LotFrontage,training.processed$log10.LotArea, use = "complete.obs")
## [1] 0.6540587
cor(training.processed$LotFrontage,training.processed$GarageArea, use = "complete.obs")
## [1] 0.3449967
cor(training.processed$LotFrontage,training.processed$MasVnrArea, use = "complete.obs")
## [1] 0.1934581
# Note that correlation gets better when log10 transformed:
cor(log10(training.processed$LotFrontage), training.processed$log10.LotArea, use = "complete.obs")
## [1] 0.7455501
summary(lm(training.processed$Log.SalePrice ~ log10(training.processed$LotFrontage) + training.processed$log10.LotArea ))
## 
## Call:
## lm(formula = training.processed$Log.SalePrice ~ log10(training.processed$LotFrontage) + 
##     training.processed$log10.LotArea)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.51402 -0.23961 -0.04072  0.25528  1.22485 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)
## (Intercept)                            8.88631    0.20139  44.125  < 2e-16
## log10(training.processed$LotFrontage)  0.30471    0.10638   2.864  0.00425
## training.processed$log10.LotArea       0.65296    0.07572   8.623  < 2e-16
##                                          
## (Intercept)                           ***
## log10(training.processed$LotFrontage) ** 
## training.processed$log10.LotArea      ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3763 on 1198 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.1833, Adjusted R-squared:  0.1819 
## F-statistic: 134.4 on 2 and 1198 DF,  p-value: < 2.2e-16
plot(y = log10(training.processed$LotFrontage), x = training.processed$log10.LotArea, cex = 0.5, col = "navy" , pch = 19)
abline(lm(log10(training.processed$LotFrontage) ~ training.processed$log10.LotArea))

It looks like this feature pretty much collinear with the LotArea feature, especially when log10 transformed. Therefore, we will remove this feature from both training and test data sets.

training.processed <- dplyr::select(training.processed, -LotFrontage)
test.processed <- dplyr::select(test.processed, -LotFrontage)

4.6 MasVnrArea: Masonry veneer area in square feet

Looking at more closely to this feature:

plot(x = training.processed$MasVnrArea, y = training.processed$Log.SalePrice,cex = 0.5, col = "navy" , pch = 19)

Note that a large number of these values are zero, which probably reflects the houses with no masonry veener. What happens if we put these aside and look at the rest of the values:

plot(x = log10(training.processed$MasVnrArea[training.processed$MasVnrArea != 0]), y = training.processed$Log.SalePrice[training.processed$MasVnrArea != 0],cex = 0.5, col = "navy" , pch = 19)

There might be ‘some’ relationship with the response, if not too strong. It is also a little confusing since the variable MasVnrType only has these 8 values as “none”, but linked the large amount of MasVnrAre == 0 data to certain MasVnrTypes:

qplot(data = training.processed,x = log10(training.processed$MasVnrArea), y = training.processed$Log.SalePrice, col = train$MasVnrType)

It is more clear now. The “none” group indeed represent the MasVnrArea = 0, NAs however appear as missing values in both cases.

ggplot(data = training.processed,aes(x = train$MasVnrType, y = training.processed$Log.SalePrice, fill = train$MasVnrType))+geom_boxplot()

We note that the SalePrice of the missing cases in the training set is closest to the SalePrice of the houses with Stone Masonry Veener type. It is likely that the area of the veener in some of these houses were not measured, perhaps it is more difficult to do so. This is our assumption.

Based on this assumption, we will merge the NA cases with the Stone group and remove the MasVnrArea feature in both training and test data sets:

training.processed$MasVnrType.None[training.processed$NO.MasVnrType == 1] <- 1
training.processed <- dplyr::select(training.processed, -NO.MasVnrType, -MasVnrArea)

test.processed$MasVnrType.None[test.processed$NO.MasVnrType == 1] <- 1
test.processed <- dplyr::select(test.processed, -NO.MasVnrType, -MasVnrArea)

4.7 Missing data remaining in the test set

We still need to deal with few missing data points we have in the test set. Since we need to make SalePrice predictions for each of these cases, we need to assess how we can impute these values at this point:

apply(is.na(test.processed), 2, sum)[apply(is.na(test.processed), 2, sum) != 0]
##   BsmtFinSF1   BsmtFinSF2    BsmtUnfSF  TotalBsmtSF BsmtFullBath 
##            1            1            1            1            2 
## BsmtHalfBath   GarageCars   GarageArea 
##            2            1            1

Some features regarding Basement and Garage have missing values:

test.processed[which(apply(is.na(test.processed), 1, sum) != 0),which(apply(is.na(test.processed), 2, sum) != 0 )]

Fortunately, the missing values are only restricted to 3 house records. It is likely these unusual houses are either lacking basement or garage. Since all these features are numeric, we will simply place zero instead of NA’s for these values. This is of course based on our assumption.

for (i in 1:ncol(test.processed)) {
        if(any(is.na(test.processed[,i]))){
                test.processed[is.na(test.processed[,i]),i] <- 0
        }
}

This completes missing value imputation as no missing values left in either training and test data sets.

5 Predictive modeling using the training set

5.1 Using shrinkage methods (regularization)

It would be interesting to perform a linear model selection and see how well it might perform with the new training set we obtained.

5.1.1 Fitting a lasso regression model to data

library(glmnet)

# We create a matrix of predictors
x.training = model.matrix(Log.SalePrice ~ . -1, data = training.processed)
x.testing = model.matrix( ~ . -1, data = test.processed)

y = training.processed$Log.SalePrice

# Fit lasso regression to training data
fit.lasso <- glmnet(x.training,y, family = "gaussian")
plot(fit.lasso, xvar = "lambda", label = TRUE)

plot(fit.lasso, xvar = "dev", label = TRUE)

Lasso is a nice way of reducing the features while still trying to explain the variance and therefore trying to win the bias-variance trade-off. In this case we notice that only 3 features are able to explain 60% of the deviance, while adding 13 features is able to explain 80% of it. It is more feasible to use the restricted model with the few features and their shrunken coefficients.

5.1.2 Choosing the optimal lambda by using cross-validation

It is best to choose the optimal lambda using cross-validation, with the aim of minimizing mean squared error::

#cv.lasso will be the list containing the coefficients and information of our optimal model fit using cross-validation (10 fold by default)
cv.lasso <- cv.glmnet(x.training,y, family = "gaussian")
plot(cv.lasso)

Notice that what we are doing here is actually computationally intense.

  1. We get very fine grids of lambda values.
  2. Using each lambda value, we perform 10-fold cross validation:
  • We fit a model using the Sum of Squares convex optimization and obtain cross-validated error (average of 10 model fits)
  1. We continue steps 1 and 2 for all fine grids of lambda values in the range.

The two vertical lines are produced in the plot, the left one marks the model with min error and the middle one shows a little more restricted model 1 standard error away from the minimum error.

When choosing the optimal model, we might prefer to get the more restricted model that is indicated by the second vertical line, which is default in glmnet package, due to the 1 standard error conversion.

The final (optimal) model coefficients we obtained by cross-validation and the non-zero features remaining:

coef(cv.lasso)
## 283 x 1 sparse Matrix of class "dgCMatrix"
##                                       1
## (Intercept)                5.309391e+00
## MSSubClass                 .           
## OverallQual                9.324533e-02
## OverallCond                1.040886e-02
## YearBuilt                  1.182768e-03
## YearRemodAdd               1.355597e-03
## BsmtFinSF1                 1.629126e-05
## BsmtFinSF2                 .           
## BsmtUnfSF                  .           
## TotalBsmtSF                5.537271e-05
## X1stFlrSF                  .           
## X2ndFlrSF                  .           
## LowQualFinSF               .           
## GrLivArea                  1.803022e-04
## BsmtFullBath               1.643841e-02
## BsmtHalfBath               .           
## FullBath                   .           
## HalfBath                   .           
## BedroomAbvGr               .           
## KitchenAbvGr               .           
## TotRmsAbvGrd               .           
## Fireplaces                 9.766781e-03
## GarageYrBlt                .           
## GarageCars                 6.852253e-02
## GarageArea                 .           
## WoodDeckSF                 1.644715e-05
## OpenPorchSF                .           
## EnclosedPorch              .           
## X3SsnPorch                 .           
## ScreenPorch                .           
## PoolArea                   .           
## MiscVal                    .           
## MoSold                     .           
## YrSold                     .           
## log10.LotArea              1.533081e-01
## MSZoning.C..all.          -1.476063e-01
## MSZoning.FV                .           
## MSZoning.RH                .           
## MSZoning.RL                .           
## MSZoning.RM               -4.101732e-02
## Street.Grvl                .           
## Street.Pave                .           
## Alley.Grvl                 .           
## Alley.NoAlley              .           
## Alley.Pave                 .           
## LotShape.IR1               .           
## LotShape.IR2               .           
## LotShape.IR3               .           
## LotShape.Reg               .           
## LandContour.Bnk            .           
## LandContour.HLS            .           
## LandContour.Low            .           
## LandContour.Lvl            .           
## Utilities.AllPub           .           
## LotConfig.Corner           .           
## LotConfig.CulDSac          .           
## LotConfig.FR2              .           
## LotConfig.FR3              .           
## LotConfig.Inside           .           
## LandSlope.Gtl              .           
## LandSlope.Mod              .           
## LandSlope.Sev              .           
## Neighborhood.Blmngtn       .           
## Neighborhood.Blueste       .           
## Neighborhood.BrDale        .           
## Neighborhood.BrkSide       .           
## Neighborhood.ClearCr       .           
## Neighborhood.CollgCr       .           
## Neighborhood.Crawfor       1.844449e-02
## Neighborhood.Edwards       .           
## Neighborhood.Gilbert       .           
## Neighborhood.IDOTRR        .           
## Neighborhood.MeadowV       .           
## Neighborhood.Mitchel       .           
## Neighborhood.NAmes         .           
## Neighborhood.NoRidge       .           
## Neighborhood.NPkVill       .           
## Neighborhood.NridgHt       5.162883e-04
## Neighborhood.NWAmes        .           
## Neighborhood.OldTown       .           
## Neighborhood.Sawyer        .           
## Neighborhood.SawyerW       .           
## Neighborhood.Somerst       .           
## Neighborhood.StoneBr       .           
## Neighborhood.SWISU         .           
## Neighborhood.Timber        .           
## Neighborhood.Veenker       .           
## Condition1.Artery          .           
## Condition1.Feedr           .           
## Condition1.Norm            4.167761e-03
## Condition1.PosA            .           
## Condition1.PosN            .           
## Condition1.RRAe            .           
## Condition1.RRAn            .           
## Condition1.RRNe            .           
## Condition1.RRNn            .           
## Condition2.Artery          .           
## Condition2.Feedr           .           
## Condition2.Norm            .           
## Condition2.PosA            .           
## Condition2.PosN           -5.019406e-02
## BldgType.1Fam              .           
## BldgType.2fmCon            .           
## BldgType.Duplex            .           
## BldgType.Twnhs             .           
## BldgType.TwnhsE            .           
## HouseStyle.1.5Fin          .           
## HouseStyle.1.5Unf          .           
## HouseStyle.1Story          .           
## HouseStyle.2.5Unf          .           
## HouseStyle.2Story          .           
## HouseStyle.SFoyer          .           
## HouseStyle.SLvl            .           
## RoofStyle.Flat             .           
## RoofStyle.Gable            .           
## RoofStyle.Gambrel          .           
## RoofStyle.Hip              .           
## RoofStyle.Mansard          .           
## RoofStyle.Shed             .           
## RoofMatl.CompShg           .           
## RoofMatl.Tar.Grv           .           
## RoofMatl.WdShake           .           
## RoofMatl.WdShngl           .           
## Exterior1st.AsbShng        .           
## Exterior1st.AsphShn        .           
## Exterior1st.BrkComm        .           
## Exterior1st.BrkFace        .           
## Exterior1st.CBlock         .           
## Exterior1st.CemntBd        .           
## Exterior1st.HdBoard        .           
## Exterior1st.MetalSd        .           
## Exterior1st.Plywood        .           
## Exterior1st.Stucco         .           
## Exterior1st.VinylSd        .           
## Exterior1st.Wd.Sdng        .           
## Exterior1st.WdShing        .           
## Exterior2nd.AsbShng        .           
## Exterior2nd.AsphShn        .           
## Exterior2nd.Brk.Cmn        .           
## Exterior2nd.BrkFace        .           
## Exterior2nd.CBlock         .           
## Exterior2nd.CmentBd        .           
## Exterior2nd.HdBoard        .           
## Exterior2nd.ImStucc        .           
## Exterior2nd.MetalSd        .           
## Exterior2nd.Plywood        .           
## Exterior2nd.Stone          .           
## Exterior2nd.Stucco         .           
## Exterior2nd.VinylSd        .           
## Exterior2nd.Wd.Sdng        .           
## Exterior2nd.Wd.Shng        .           
## MasVnrType.BrkCmn          .           
## MasVnrType.BrkFace         .           
## MasVnrType.None            .           
## MasVnrType.Stone           .           
## ExterQual.Ex               .           
## ExterQual.Fa               .           
## ExterQual.Gd               .           
## ExterQual.TA               .           
## ExterCond.Ex               .           
## ExterCond.Fa               .           
## ExterCond.Gd               .           
## ExterCond.Po               .           
## ExterCond.TA               .           
## Foundation.BrkTil          .           
## Foundation.CBlock          .           
## Foundation.PConc           .           
## Foundation.Slab            .           
## Foundation.Stone           .           
## Foundation.Wood            .           
## NO.BsmtQual                .           
## BsmtQual.Ex                2.446339e-02
## BsmtQual.Fa                .           
## BsmtQual.Gd                .           
## BsmtQual.TA                .           
## NO.BsmtCond                .           
## BsmtCond.Fa                .           
## BsmtCond.Gd                .           
## BsmtCond.Po                .           
## BsmtCond.TA                .           
## NO.BsmtExposure            .           
## BsmtExposure.Av            .           
## BsmtExposure.Gd            9.061738e-04
## BsmtExposure.Mn            .           
## BsmtExposure.No            .           
## NO.BsmtFinType1            .           
## BsmtFinType1.ALQ           .           
## BsmtFinType1.BLQ           .           
## BsmtFinType1.GLQ           8.339921e-03
## BsmtFinType1.LwQ           .           
## BsmtFinType1.Rec           .           
## BsmtFinType1.Unf          -6.834311e-03
## NO.BsmtFinType2            .           
## BsmtFinType2.ALQ           .           
## BsmtFinType2.BLQ           .           
## BsmtFinType2.GLQ           .           
## BsmtFinType2.LwQ           .           
## BsmtFinType2.Rec           .           
## BsmtFinType2.Unf           .           
## Heating.GasA               .           
## Heating.GasW               .           
## Heating.Grav               .           
## Heating.Wall               .           
## HeatingQC.Ex               1.334368e-02
## HeatingQC.Fa               .           
## HeatingQC.Gd               .           
## HeatingQC.Po               .           
## HeatingQC.TA               .           
## CentralAir.N              -6.204003e-02
## CentralAir.Y               .           
## Electrical.FuseA           .           
## Electrical.FuseF           .           
## Electrical.FuseP           .           
## Electrical.SBrkr           .           
## KitchenQual.Ex             4.700208e-03
## KitchenQual.Fa             .           
## KitchenQual.Gd             .           
## KitchenQual.TA            -1.387015e-02
## Functional.Maj1            .           
## Functional.Maj2            .           
## Functional.Min1            .           
## Functional.Min2            .           
## Functional.Mod             .           
## Functional.Sev             .           
## Functional.Typ             .           
## FireplaceQu.Ex             .           
## FireplaceQu.Fa             .           
## FireplaceQu.Gd             .           
## FireplaceQu.NoFireplace   -2.957736e-02
## FireplaceQu.Po             .           
## FireplaceQu.TA             .           
## NO.GarageType              .           
## GarageType.2Types          .           
## GarageType.Attchd          6.147442e-03
## GarageType.Basment         .           
## GarageType.BuiltIn         .           
## GarageType.CarPort         .           
## GarageType.Detchd          .           
## NO.GarageFinish            .           
## GarageFinish.Fin           .           
## GarageFinish.RFn           .           
## GarageFinish.Unf           .           
## NO.GarageQual              .           
## GarageQual.Fa              .           
## GarageQual.Gd              .           
## GarageQual.Po              .           
## GarageQual.TA              .           
## NO.GarageCond              .           
## GarageCond.Ex              .           
## GarageCond.Fa              .           
## GarageCond.Gd              .           
## GarageCond.Po              .           
## GarageCond.TA              .           
## PavedDrive.N              -5.958627e-03
## PavedDrive.P               .           
## PavedDrive.Y               .           
## PoolQC.Ex                  .           
## PoolQC.Gd                  .           
## PoolQC.NoPool              .           
## Fence.GdPrv                .           
## Fence.GdWo                 .           
## Fence.MnPrv                .           
## Fence.MnWw                 .           
## Fence.NoFence              .           
## MiscFeature.Gar2           .           
## MiscFeature.NoMiscFeature  .           
## MiscFeature.Othr           .           
## MiscFeature.Shed           .           
## SaleType.COD               .           
## SaleType.Con               .           
## SaleType.ConLD             .           
## SaleType.ConLI             .           
## SaleType.ConLw             .           
## SaleType.CWD               .           
## SaleType.New               .           
## SaleType.Oth               .           
## SaleType.WD                .           
## SaleCondition.Abnorml     -7.582707e-03
## SaleCondition.AdjLand      .           
## SaleCondition.Alloca       .           
## SaleCondition.Family       .           
## SaleCondition.Normal       .           
## SaleCondition.Partial      .

Note that our final model includes ~30 features with non-zero coefficients. We therefore significantly reduced the size of the data set by selecting features through lasso regression. The resulting model explains more than 80% of the deviance in the outcome, which is considerable good for the training set.

We also note that some of the dummy variables remained in the model, while majority of them removed. It seems to be good idea to keep all levels as seperate dummy variables until the regularization, since we don’t know which ones may be important in the final model context.

5.1.3 Making predictions using the best lasso model fit

Let’s first start with the training set:

pred <- predict(cv.lasso, newx = x.training)

#Calculating the root mean squared error
rmse.training = sqrt(mean((training.processed$Log.SalePrice - pred)^2))
rmse.training
## [1] 0.1479088
#Pearson correlation:
cor.training = cor(training.processed$Log.SalePrice,pred)

plot(x = pred, y = training.processed$Log.SalePrice, cex = 0.5, col = "navy", pch = 19)
text(x = max(pred), y = min(training.processed$Log.SalePrice)+0.3,
     labels = paste0("RMSE: ",rmse.training))

Therefore, the RMSE we obtained from the training set is ~15%.

Next we will try to use test set to estimate the house prices using our model.